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Abstract. Consider an ideal I C R ~ C[xi, . . . , Xn] defining a complex affine variety 
X C C". We describe the components associated to / by means of numerical primary 
decomposition (NPD). 

The method is based on the construction of deflation ideal l'^'') that defines the deflated 
variety X^'^^ in a complex space of higher dimension. For every embedded component 
there exists d and an isolated component Y^'^'> of 7^''-' projecting onto Y. In turn, Y^'^'> can 
be discovered by existing methods for prime decomposition, in particular, the numerical 
irreducible decomposition, applied to X^'^\ 

The concept of NPD gives a full description of the scheme Spec(i?//) by representing 
each component with a witness set. We propose an algorithm to produce a collection of 
witness sets that contains a NPD and that can be used to solve the ideal membership 
problem for I. 

1. Introduction 

Throughout the paper we use the following notation. Let / be an ideal in the polynomial 
ring R = C[x] = C[a;i, . . . , generated by polynomials fi,...,fN £ R- The ideal 
/ defines the a variety X = V{I), which, set-theoretically, is the set of points in C" 
annihilated by all polynomials in ideal /. 

The basic description of the affine scheme Spec(-R//) is given by the set of associated 
prime ideals Ass(/) consisting of all ideals that appear as annihilators of elements of the 
i?-module R/I. 

There exist symbohc algorithms (see, e.g., [SI El El El [13 [IS] ) to compute an irredundant 
primary decomposition of /, which is by definition a decomposition 

(1.1) I = Jin---nJr 

such that all ideals Jj are primary and their radicals li = \fTi G Ass(/) are pairwise 
distinct. 

The set of subvarieties defined by the prime ideals 

VAss(/) = {V{P) I P G Ass(/)} 

is called the components associated to /. A component is said to be isolated if it is maximal 
with respect to inclusion and embedded if it is not. The existing methods of numerical 
algebraic geometry can "see" only the isolated components. 

The concept of numerical primary decomposition (NPD) that we propose describes all 
components in terms of the (generalized) witness sets (see Definition 14. ip . Such a witness 
set for a component Y G VAss(J) describes Y completely: one can sample the component, 
determine its degree, determine whether a given point belongs to Y, etc. Moreover, with a 
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NPD in hand one can solve the ideal membership problem, i.e., given a polynomial f & R 
decide if / G /. This is discussed in Subsection 14.21 

The idea of our method lies in the construction of the deflation ideal I^^^ (see Definition 
12.31) that defines the deflated variety X^'^^ = V{I^'^^) of order c? in a higher-dimensional 
space. The latter is a stratified vector bundle over X that comes with the natural projec- 
tion lid : X*^'^'' — s> X onto the base. This construction is related to a deflation technique 
for 0-dimensional isolated components [121 1131 E] ■ 

The main theoretical result, Theorem 13.81 shows that an embedded component Y G 
VAss(/) becomes visible for some deflation order d, i.e., there is an isolated component 
Z e VAss(/('^)) such that iraiZ) = Y. 

Based on this theorem, we outline a straightforward algorithm. Algorithm 13. IH to 
compute all visible components up to the given order. The advantage of this algorithm 
over other known techniques is that the problem of finding all components is reduced to 
that of finding only isolated components without performing any saturation steps. 

Algorithm 15.31 specializes this general algorithm to a numerical method for prime de- 
composition, namely, the numerical irreducible decomposition (NID) |17j . 

The area of numerical algebraic geometry comprises novel approaches to computational 
algebraic geometry based on the numerical polynomial homotopy continuation methods. 
The recent book fTEl may serve as a good introduction to the area. 

While these methods involve computations that are approximate, a typical output of 
the algorithms consists of not only approximations of exact solutions, which are algebraic 
numbers, but also exact discrete information about the input. For instance, the main 
concept introduced in this paper, a numerical primary decomposition of an ideal provides 
exact data such as the dimensions and the degrees of the associated components. Another 
example is a solution to the ideal membership problem given via NPD. Although we do 
not provide a certification procedure for this method, such can be developed in theory, 
provided that results of our numerical computation can be refined to an arbitrary precision. 

One advantage of hybrid numerical homotopy continuation techniques over purely sym- 
bolic methods, such as Grobner bases, is that the former is easily parallelizable and the 
algorithms for computing the latter are intrinsically serial. 

Due to a very small amount of data being stored on or transferred between the com- 
putational nodes and embarrassing parallelism, one can easily achieve linear speedups for 
numerical homotopy continuation on a parallel computing system with any architecture. 
This comes in contrast with the high interdependency of tasks performed in Grobner bases 
computation and the phenomena of intermediate expression swell. (Although there is no 
doubt that limited speedups are possible in parallel Grobner computation, the claims of 
good scalability are substantiated with experiments on preselected classes of problems and 
methods based on non-optimal serial algorithms. See, e.g., [1] for critique of such claims.) 

We live in the age when the clock speed of processors has stopped growing fast and the 
computational capacity of computers increases mainly through building either multicore 
or distributed systems. This tendency implies both great present and even better future 
for the algorithms of numerical algebraic geometry. 



NUMERICAL PRIMARY DECOMPOSITION 



3 



The paper is structured as follows. After the introduction (Section 1), we make the 
main definitions of the paper - those of deflated ideal and variety - in Section 2. Next, 
we overview the dual space approach to looking at the structure of a polynomial ideal; 
Section 3 culminates in the proof of Theorem I3.8[ the main theoretical result of the paper. 
Section 4 proposes a new numerical representation of an ideal called numerical primary 
decomposition (NPD) via (generalized) witness sets. In Section 5, we give a skeleton of 
an NPD algorithm and examples. We discuss the future of growing meat on the skeleton 
in conclusion, Section 6. 



2. Deflation 

A variant of the matrix defined below appears in the deflation method for regulariz- 
ing singular isolated solutions of a polynomial system described in [13], as well as the 
computation of the multiplicity of an isolated point [1]. 

Definition 2.1. The deflation matrix of order d of an ideal I generated by polynomials 
fi,...,fN^Risa matrix A^f^ with 

• entries in R; 

• rows indexed by x°'fj, where \a\ < d and j = 1, 2, . . . , N; 

• columns indexed by partial differential operators = — 3^— -a;^, where \/3\ < d; 
and with the entry at row x'^fj and column set to be 

(2.1) . [xV)) = 

If a point X is an isolated solution, in other words, {x} is 0-dimensional irreducible 
component, then its multiplicity - i.e., the vector space dimension of the ring R/I localized 
at X - equals corank A^f^ (x) for a large enough d. This fact is shown, for example, in [1] 
where a method for computing the multiplicity using A^f'\x.) is introduced. 

Remark 2.2. A variant of the deflation matrix - the column of A^f^ labeled with is 
dropped - is used in the (higher-order) deflation in [12] to construct an augmented system 
of equations for which the multiplicity of a given isolation solution x of the original system 
drops. For computational purposes there is no need to keep the aforementioned column, 
however, this paper's definition of deflation matrix makes our argument in Section [3] more 
compact. 

Definition 2.3. Let I = (/i,...,/Ar) C R and let a = (a^), \j3\ < d, be a vector of 
indeterminates. 

The ideal generated by fi, . . . , fj^j and the entries of the vector A^^^a'^ in the ring C[x, a] 
is called the deflation ideal of I of order d and denoted by I^'^^ . 
The deflated variety of order d is defined as 

where B{n, d) = dimC[a3, a] = n + {^~^'^~^) ■ 
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Given an ideal /, the variety X^'^^ is well defined in view of the following proposition, 
the proof of which is self-contained. 

Lemma 2.4. For every g E I , the deflated variety Y^'^'^ of its hypersurface Y = V{g) 
contains X^'^'^ = V{I^'^^). 

Proof. Let Qa = Yl ^^pd^ ^ C[a, d] be the linear differential operator corresponding to 
the vector a = (a/3) of indeterminates. It is enough to show that for every point x G X 
and every vector a in the kernel of A^^\x) the expression Qag vanishes at x. This would 
imply that the fiber of Y^'^^ Y over x contains the fiber of X^^^ X over x, thus 
proving the statement. 

Let us write g = J2ai<i<N ^h^x" fi, where Ci^a £ C. Given a point x G C" we can 
rewrite this as two sums: 

9= E <«(^-x)"/.+ <"^"/- 

\a\>d,l<i<N \a\<d,l<i<N 

Now, apply the operator Qa, to g. The first sum above vanishes at x, since the order 
of vanishing of every summand is at least d + 1. The second sum vanishes, since a G 
kevAfi^). □ 

Proposition 2.5. If h C h C R, then K(/f ^) C 1/(/f^) for all d. 

In particular, the deflation variety X^'^^ does not depend on the set of generators fi, . . . , f^ 
of I chosen to construct Af"^ . 

Proof. This follows from Lemma [2741 since for a fixed set of generators /i, . . . , /jv of /, by 
definition, 

V (/(^)) = vii) n V ((/i)(^)) n---nv {{f^Y"^) . 

□ 

Moreover, a stronger Proposition 12.71 holds; its proof requires a global argument less 
transparent than the local argument of Lemma 12.41 and 12. 5[ 

Lemma 2.6. Let Qa = X]|/3|<d '^/3^^ ^'^^ XiQa is a formal derivative of Qa with respect 
to di for 1 < i < n. 

Then x"Qa ■ fj G I^'^^ for every j and any a, where the deflation ideal I^'^^ is defined 
using fi,...,fN- 

Proof. In case |a| > d, the expression x^Qa is a constant in d, hence, x^'Qa ■ fj £ I^'^\ 
since fj G I^'^\ 

Assume \a\ < d. By generalized Leibnitz rule 

(2.2) = E 7^^^^-/- 

7 

The expression is zero if not 7 = a or I7I < Assuming that the claim is 

established for degrees lower than it follows that x?Qa ■ fj £ -^^'^^ for I7I < Since 
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4)(x)= 

XIX2X3 



the left hand side of ( 12.2p belongs to I^'^^ by definition, it follows that x'^Qa ■ fj, the 
summand corresponding to 7 = a, belongs to I^^^ as well. □ 

Proposition 2.7. The deflation ideal J*^'^-* does not depend on the set of generators I 
chosen to construct it. 

Proof. By Lemma [2.6^ we have Qa ■ (x^fj) G therefore, Qa ■ h G I^'^^ for every h G I, 
since /i is a linear combination of terms x'^fj and the action of Qa is linear. 

The conclusion of the proposition follows immediately. □ 

Example 2.8. Let / = (xf, xiX2X-i) C C[xi, X2, X3]. The radical y/l = (xi), hence, the 
only isolated component is V{xi). 

Compute the first order {d = 1) deflation: multiply the deflation matrix 

id di 82 ^3 
^ xf 2x1 ^ 

X1X2X3 X2X3 X1X3 XiX2_ 

by the column vector (ao, cti, 02, a^)^. The deflation ideal then is 

J^^^ = {xl, X1X2X3, aiXi, 01X2X3 + 02X1X3 + 03X1X2). 

Observe that 

v/w = (xi, 01X2X3) = (xi, oi) n (xi, X2) n (xi, X3), 

therefore, projecting isolated components of X'^^^ to X gives two more (embedded) com- 
ponents: y(xi,X2) and l^(xi,X3). 

Remark 2.9. Denote by tt^ : X^'^'^ X the natural projection induced by the projection 
71 d '■ C^^^''^^ C", which maps (x, a) t-^ x. 

The deflated variety X^'^^ is a stratified vector bundle over the variety X with the pro- 
jection Tid onto the base. The strata of X over which X^'^'' is locally trivial are constructible 
algebraic subsets of X. 

Example 2.10. Let L be a plane in C" of dimension n — k defined by the vanishing of 
/ = (xjj for some 1 < Ji < ■ ■ ■ < Ja; < 
Then it is not hard to establish that 

= / + {{ap ■.\(3\<d and f3j^ ^ for all i}). 

Therefore, L^'^^ is a plane in C^'^"''^^ of dimension n — /c + (" ^^'^ ^) . 

Every plane in can be brought to the above form with an affine change of coordinates. 
This change would result in a linear transformation on the derivations d^; hence, a deflated 
variety of order d of a plane of codimension k is always a plane of dimension n — k + 

3. Visible components 

In this section we describe briefly the dual space (a.k.a. inverse system) approach to 
the local description of a point scheme. While from the computational point of view this 
approach has not been exploited yet beyond the 0-dimensional case, it comes handy in 
the proof of the main theorem of the section and the paper. 
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3.1. Dual spaces. 

Definition 3.1. For X G C", let Al^: R^C be a lin ear functional defined by 

A^(/) = (^"-/)(x) = ^(x), feR. 

The dual space -Dx[-^] of ideal I at a point x is the subspace of the C-span of functional 
consisting of the ones that annihilate all polynomials in I. 

The dual space defined as above has been introduced under other names in literature; 
notably it is the local inverse system, the term going back to Macaulay [Hj. 
The dual space -Dx[-^] has a filtration 

[/]cdW[/]cD(^)[/]c... 

where D^^ [I] is the set of functionals of order at most d. For computational purposes, it is 
convenient to think of D^^ [I] as ker Af"^ (x) : indeed, given a vector a e ker A^'^ (x) , one can 
view it naturally as the functional Yl\^\<d^P^x- The natural isomorphism ker A^'^^(x) ~ 
D^\l] could be shown then by inspection of the definition of the defiation matrix. In 

particular, for the case of an isolated point x, the whole dual space D^[I] ~ ker A^'^-' (x) 
for d ^ 0, which is shown, for example, in [U Theorem 1]. 

Let Ci be an ra- vector with one in the z-th position and zeros in the rest. On the space 
D^[0] of all differential functionals at a point x we define 

(1) operators of integration 6i : -Dx'^''[0] D^^^^O], i = induced by the 
multiplication by di in the Definition 13.11 i.e.: 

6,{Ai){f) = {d,d^-f), feR, 

in other words, 

(2) operators of differentiation Xi '■ D^^^\o] — > Dx'^^[0], i = 1, . . . ,n, induced by the 
multiplication by (xj — Xj) in the Definition 13. H i.e.: 

X.(A^)(/) = ((x,-x,)9^-/) 

= {{d^{x.-^,)+(3,d^-'^)-f), 

therefore, after specializing a; = x, 

X.(A^) = AA^^^ 

omce A^ = S^A°, for a fixed X we may look at the space of linear functionals -Dx[0] 
as C[(5] = C[6i, . . . , 6n\ with differentiation operators Xi acting as formal derivations 
Note that -Dx[-^] C C[(5i, . . . , 5„] is stable under this action of Xi for every ideal I <Z R. 

In fact, by (naturally) extending the space of linear functionals to the formal power 
series C[[(5]] one may obtain the following fact (see, e.g., [151 Proposition 2.6]). 



Proposition 3.2. The ideals of R are in one-to-one correspondence with the vector sub- 
spaces ci/C[[5]] stable under differentiation and closed in {d)-adic topology. 
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For our purposes, given an ideal J, we need just a local statement regarding (_R//)x = 
{R/I)m^, the quotient local ring at a point x, where rrix = (xi — xi, . . . , — x„) C R. 

One may show that the subspaces of C[5] stable under the differentiation are in one- 
to-one correspondence with the ideals of the local ring R^. The following lemma proves 
this statement in one direction. 

Lemma 3.3. The image of f ^ R is zero 

(1) m {R/I + mf zjf Q- f = for allQ e oi^^ [I]; 

(2) zn (i?//)x zff Q-f = OforallQe D^[I] . 

Proof. Statement (2) follows from (1). 

To prove (1), notice that D:x.[I + m'^^] is exactly Dx\l]- Indeed, every functional of the 
latter annihilates not only /, but also all elements of m^+^; on the other hand, a functional 
of order larger than d does not kill the entire m^^^. □ 

Remark 3.4. In case of an isolated point x G X, it follows from the lemma that the 
dimensions of the C-spaces -Dx[-^] and (i?//)x are equal. Stronger statements are available: 
for example, see fi3{ Theorem 3.1] and [151 Theorem 3.2]. The latter provides a recipe for 
constructing the iTix-primary component of the ideal / from a basis of -Dx[-^]- 

3.2. A visible deflation of a component. Recall the deflated variety X^'^^ = V"(/'^'^^) C 
(j^s(n,d) g^j^j ^Yi^Q projection vr^ : X^'^^ — > X in the Definition 12.31 

Definition 3.5. We define the deflation of order d of a component Y G VAss(/) as 
y('^) = TT^^Y° C X^'^\ where Y° is the subset o/ generic points, that is all smooth points 
that do not belong to other components that do not contain Y . 

Proposition 3.6. A deflation of a component is an irreducible subvariety. 

Proof. As was mentioned in Remark 12.91 X^^^ is a stratified bundle. The generic locus 
Y° is a stratum in the stratification; indeed, the bundle is locally trivial at every generic 
point. Since t^^^Y" is an open subset of Y^^\ the latter is irreducible. □ 

Definition 3.7. We say that Y G VAss(/) is visible at order d, if Y^'^'' is an isolated 
component of X^'^^ . 

Theorem 3.8. Every component is visible at some order. 

Proof. Note that an isolated component Y is visible for any d > 0, since Y^'^'' = tx^^Y° 
can not be contained in any other component in VAss(/*^''''). 

Let / = flzeVAssC/) '^z be an irredundant primary decomposition of /, where V{Jz) = Z 
for all Z. Fix an embedded component Y with a generic point y G F° and consider 
^>Y = nyczGVAss{/) "^z- Then />y ^ /, moreover, the local quotient ring (_R//>y)y is a 
proper quotient of {R/I)y. Hence, according to Lemma 13.31 there is a proper inclusion 
Z?y[/]2^y[/>y]. 

For any Z G VAss(/) containing Y the fiber of its deflation tt^ ^(y) fl Z^'^^ is contained 
in Z}y'^'*[J>y], since only the fibers over the generic locus Z° define Z'^'^\ Now, there exists 
d suich. tti3jlj 
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that implies immediately that Z^"^^ does not contain Y^'^\ □ 

Remark 3.9. If y is visible at order d then it is visible at any order d' > d. Indeed, the 
component 

where Tid'^d '■ X^'^ ^ X^^^ is the natural projection, has to be isolated if Y^'^^ is isolated. 

Since there is a finite number of components there is an order, at which all components 
are visible. 

Example 3.10. Consider two ideals of C[xi, 0:2, xs]: 

/ = {xl,XiX2X3) 

J — (^]^, X\X2X'^ . 

The first order deflation of / has been computed in Example 12.81 For J we have 

id di 82 83 



X 



= Xixlx3 



x{ 2x1 

X\X2X3 X2X3 2^X^X2X3 X\X2 
X1X2X3 \_XiX2x\ X2X3 Xixl 2X1X2X3 

The entries of A^j\ao, ai, a2, 03)"^ together with the original generators of J generate 



xf, X1X2X3, xiX2x|,aiXi, 



01X3X3 + 202X1X2X3 + O3X1X2, 

OiX2x| + 02Xix| + 2O3X1X2X3 ). 

It is easy to see that 

a/jw = a/jw = (x^, 01X2X3) = (xi, oi) n (xi, X2) n (xi, xs), 

therefore, we need a higher order deflation to distinguish J from I. 

One may check that the second order deflation uncovers another embedded component 
in VAss( J), the origin, which is not a component in VAss(J). However, still V J(^) = VT^, 
for / the origin is a pseudo-component! 

It is not until the third deflation that we see the difference and here is why. The 
difference may be seen by comparing the dual spaces at the origin: Do[I] and -Do[<^]- First 
of all, Do[I] C Do[J], since J d I . This inclusion is proper as 
functional of order 3. On the other hand, D^\l\ = D^\j]. 

Here we would like to propose the most general algorithm for primary decomposition, 
which take the order of deflation as a parameter and returns all components for (i ^ 0. 

Algorithm 3.11. A/" = VisibleComponents(/, c/) 

Require: I, ideal of R = K[x] where K is a field of characteristic 0; d > 0. 
Ensure: C, the set components visible at order d. 

= 71(1 {^isolated components of I^'^^^. 
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The correctness of the algorithm is assured by Theorem 13. 8t although we use K = C, 
the argument holds for an arbitrary algeraically closed field of characteristic 0. 

In fact, by passing to the algebraic closure, one can make the same conclusion for an 
arbitrary K of characteristic 0. The only serious issue in this case is that a component 
that is irreducible over K may become reducible over its algebraic closure. However, the 
deflation procedure preserves irreducibility and the property of algebraic closedness is 
exploited only locally in the proof of Theorem 13.81 

While any routine for prime decomposition can be used to compute the isolated com- 
ponents of the deflated variety, in what follows we concentrate on a numerical approach, 
which is applicable only for IK = C. 

4. Witness sets and numerical primary decomposition 

4.1. Generalized and classical witness sets. All numerical algorithms based on ho- 
motopy continuation boil down to the computation of approximations to points of a 
0-dimensional variety. That is why for every component Y we need to invent a presenta- 
tion that would consist of a finite number of points and, perhaps, some additional (finite) 
information. 

Definition 4.1. A witness set W = Wy of a component Y G VAss(/) is a triple 
{d,L,w) = {dY,LY,WY) consisting of 

(1) an order d, such that Y is visible at order d; 

(2) a generic {codimY^'^^) -plane L C C^*^"'''-'; 

(3) the (finite) set of witness points w = Y^^^ fl L; 

All items can be presented with finite data: in particular, L can be represented by a 
linear basis. We do not include as elements of the witness set generators of /, we assume 
that those are fixed and available. 

We also assume that there is a procedure Hl,l' that for another generic (codimF^'^))- 
plane L' C C^*-"''^'' takes the witness points wy as input and produces a new set of witness 
points w'y forming a witness set (c?, L',Wy). In numerical algebraic geometry such a 
procedure is provided by a sufficiently randomized homotopy continuation that deforms 
L into L' without encountering an intermediate plane that is singular with respect to Y^'^^ 
and a numerical routine that tracks the paths starting at the witness points wy. 

We would like to remark that, in principle, it is enough to store only one witness 
point, since the rest can be obtained due to the action of the monodromy group, which 
is transitive on wy. In practice, this can be done by following a random homotopy cycle 
T-Cl,l a finite number of times. However, it is, of course, more practical to store the whole 
set Wy. 

Remark 4.2. An isolated component is visible at order d = 0. In this case, the plane L 
and the set w give what we would call a classical witness set, a concept which is used, 
for example, throughout [18] where a generalization of it is made (page 237) in relation 
to the defiation of an isolated (but multiple) component. 
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Example 4.3. The components of / = Xz{x2 + x\)) C C[xi, 0:2, xs] are Z = V^xs) 
and Y = V{x3, X2 + xf). 

The isolated component Z can be presented by the witness set using any line L that is 
not parallel to Z and such that Z Cl L ^ Y . 

To represent Y, we have to look at the first order deflation: 

= (xg, Xsix2+xl), 03X3, 

2aiXiX3 + 023^3 + a'3ix2 + xl)) 

C C[xi,X2,X3,ao, 01,02,03] 
X« = Vixs, 03(^2 + C 
The deflated variety X^^^ is 5-dimensional; we take the following 5 equations for L: 

X2 = -3xi + 2; 
O/3 = Xi — 85 

Oj = CjXi + di, Ci, di E C, i = 0, 1, 2. 

The first two equations together with the second defining equation of X^^'^ give 

{xi - 3){xi - 2){xi - 1) = 0. 

The set 7ri(X'^^)nL) = {(3, —7, 0), (2, —4, 0), (1, —1, 0)} contains projections of two subsets 
of witness points, 

7ri(^i;z) = {(3,-7,0)}, 
7ri(w) = {(2, -4,0), (1,-1,0)}, 
of the witness sets of the first order Wz = (1, L, wz) and Wy = (1, L, wy), respectively. 

Remark 4.4. For a witness point (y, a) G wy, the vector a translates into a functional 
Qa £ Dy\l]. For a fixed generic y G F the set {Qa \ (y, «.) G F^*^)} equals the dual space 
D^\l]. Therefore, in practice, we can compute D^\l] from the witness set {d,L,w) of 
Y by tracking a homotopy that creates another witness set {d, Li, Wi), where the plane Li 
is a random plane such that 7rd(L) DY 3 y. If the procedure is carried out for sufficiently 
many Lj, then the functionals {Qa \ (y? «-) G Wi} span Dy'^'*[J]. 

4.2. Numerical primary decomposition and the ideal membership problem. 

Definition 4.5. A collection of witness sets is called a numerical primary decomposition 
(NPD) of I if it contains precisely one witness set for each component in VAss(/). 

NPD contains exhaustive information about the ideal / and, in particular, the scheme 
Spec(i?/J) due to the possibility of solving the ideal membership problem. This is so, since 
for every Y the projection 7ifiy{z) of a witness point z G wy gives a generic point of Y 
and in the view of the following: 

Theorem 4.6. A polynomial g & R is contained in the ideal I iff for all Y G VAss(/) and 
every (any) generic pointy G Y, all functionals in the dual space D^'^^^^\l] annihilate g. 
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Proof. A polynomial e / iff its image in i?// is zero or, equivalently, for all x G X = V{I) 
its image in {R/ I)y^ vanishes. It suffices to check the latter statement for one generic point 
per component for all components. □ 

It follows that, in practice, we can solve the ideal membership problem by checking the 
condition in the theorem at a finite number of points. Therefore, we have the following 
algorithm. 

Algorithm 4.7. h = lMF{gJ,Af) 

Require: I, ideal of R represented by a finite set of generators; N" , a NPD of I. 
Ensure: b = "(7 G a boolean value. 

Let d = degg and g = E|/3|<d C/J^;^; 
for all {d', L,w) e Af do 
Pick X G T^'diw); 

Compute a linear basis K 0/ ker Aj'^^ (x) ; 

if Q ' f ^ (*) /'^'^ some Q G D^^ [I] corresponding to an element of K then 

Return false; 
end if 
end for 

Return true. 

According to Theorem 14.61 / G / iff there is no x G X, for which the condition (*) holds. 
In fact, it suffices to check (*) for a set of generic points (one per component). 

Remark 4.8. In case the above algorithm is executed numerically., i.e., only approxi- 
mations of points on the components are generated, we would like to point out several 
practical issues. 

First, having the whole NPD it easy to generate other points on any given compo- 
nent and recheck the condition (*) at as many points as desired, therefore, lowering the 
probability of this algorithm returning an incorrect result due to picking a non-generic 
point. 

Second, since the approximations of generic points can be refined to an arbitrary pre- 
cision, the condition (*) can be effectively checked by "zooming in" on the exact point 
for which the computation is carried out. However, a rigorous certification procedure has 
to be developed in the future in the spirit of the alpha-test for an approximate zero of a 
univariate polynomial. 

Algorithm 14. 7l is practical only for polynomials of low degrees due to the high complexity 
of construction of the deflation matrix of order d and computing its kernel. However, we 
are confident that improvements can be made as this matrix is highly structured and the 
fact that it is enough to check the condition (*) for only one generic element in the kernel. 

5. Algorithm for numerical 

PRIMARY decomposition 

In the description of our NPD algorithm, we assume the following subroutines are at 
our disposal: 
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Algorithm 5.1. M = NID(/) 
Require: I, ideal of R. 

Ensure: Ai, the set of classical witness sets (Ly,wy) for all isolated components Y G 
VAss(J). 

There are two approaches to NID in the numerical algebraic geometry: the "top-down" 
method described in detail in [HI Chapter 15] and "equation-by-equation" method the 
philosophy of which is outlined in [T8l §16.2]. 

Algorithm 5.2. b = IsInComponent(y, VF^) 

Require: y G C"; Wz = {d, L, w) a witness set for Z G VAss(/). 
Ensure: b = G Z", a boolean value. 

Pick a generic (codiia Z ) -plane M C C" and a (dim L)-plane V C C^^"-'*^) such that 
Tia{L') = M3y. 

Use a procedure similar to the usual containment test routine [181 §15.1], i.e., track the 
points w along a generic homotopy 'Hl,l'- 

Return "y G 7id{nL,L'iw))". 

There are also two subroutines, for which finding efficient algorithms is an open problem: 

• StopCriterium((i, /, A^) implements a termination criterion that guarantees that 
all components in VAss(/) are visible at order less than d. 

• IsComponent {{d, L,w), I ,Af) is used to filter out witness sets {d,L,w) that rep- 
resent false components that appear due to singularities. As a parameter it takes 
partial NPD A/" that includes the witness sets for all components visible at order 
d of dimensions higher than the dimension of the alleged component. 

Both routines will be discussed later in this section. 

Now we are ready to outline the main algorithm of this paper. 

Algorithm 5.3. Af = NPD(/) 
Require: I, ideal of R. 

Ensure: M , the set of witness sets for all Y G VAss(/) visible at order d. 

M = %; 
repeat 

U = {{L,w) G N1D(J('^)) : 

not lslnComponent(7rrf(wi), W) for all l/T G C}; 
for all (L, w) E M in an order 
of decreasing dimY(^Li,j-) do 

if IsComponcnt(((i, L, ty), /, A/") then 

Af = J\fU{{d,L,w)}; 
end if 
end for 
d = d+l; 
until StopCriterium((i, /, A/") 
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The order of deflation sufficient to discover all components can be bound by the maxi- 
mum of the regularities of the (local) Hilbert functions at the points of X. For the latter 
a crude bound exists, which is doubly exponential in the number of variables. Obviously, 
one can not use StopCriterium based on this bound for practical purposes. 

In reality, for many nontrivial examples all embedded components are discovered by 
the deflation of order as low as 1. While finding a reasonable termination criterion is the 
matter of the future, currently it makes sense to run a truncated computation with 

StopCriterium (d, /, A/") = ''d > dmax" , 

where dmax is the maximal deflation order considered. 

As to IsComponent, first of all, we remark that elimination of fake components from 
J\f may also be done at the end of the algorithm. There is a way to do this by checking 
whether for large enough d all functionals in the D^\l] for a generic y in an alleged 
component Y "come" from the components that contain Y. We do not describe the 
procedure here as it is quite technical and not practical at the moment, since the question 
how large d should be relates to the question of finding a good StopCriterium. 

Remark 5.4. Establishing an efficient StopCriterium has a higher priority (over IsComponent), 
since Af containing additional witness sets of fake components can be used instead of a 
true NPD for many tasks. In particular, our IMP routine (Algorithm 14.71) would still 
work. 

Below is an example of a "truncated computation" , where only the first order deflation 
is computed. 

Example 5.5. Consider the cyclic 4-roots problem: 

/ = (xi + X2 + X3 + X4, X1X2 + X2X3 + X3X4 + X4X1, 
X1X2X3 + X2X3X4 + X3X4X1 + X4X1X2, XiX2X^Xi — 1) . 

The calculation of associated primes via symbolic software (we used Macaulay 2 [S]) 
gives: 

Ass(/) = { (x2 + X4, xi + X3,X3a;4 + 1), 

{X2 + X4, Xi + X3, X3X4 - 1), 

{Xi - 1, X3 + 1, X2 + 1, Xi - 1), 

(X4 - 1,X3 - 1,X2 + + 1), 

(X4 + 1, X3 + 1, X2 - 1, Xi - 1), 

(X4 + 1, X3 - 1, X2 - 1, Xi + 1), 

(X3 + X4, X2 + X4, Xi - X4, X4 + 1), 

(X3 -X4,X2 +X4,Xi +X4,X4 + 1) } 

The first two ideals correspond to the irreducible curves that are the two isolated com- 
ponents. The rest are embedded 0-dimensional components; note that the last two ideals 
are irreducible over the ground field Q, but not C. 

Over complex numbers, there are 8 embedded components that are all visible at order 1; 
The numerical computation of the irreducible components of the first deflation finds all 
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components. An excerpt from the numerical output is given below: we list the projections 
of the witness points for the components VAss(/^^''). 

»> projections of witness points for 
»> component #1 : 

[xl = -4.4882+2.0260*1, x2 = . 18509+ . 83550e-l*I , 
x3 = 4.4882-2.0260*1, x4 = -. 18509- . 83550e-l*I] 
[xl = .52885e-l-. 87608*1, x2 = -. 68654e-l-l . 1373*1 , 
x3 = -.52885e-l+. 87608*1, x4 = . 68654e-l+l . 1373*1] 
[xl = -.15083+. 49191*1, x2 = .56975+1.8582*1, 
x3 = .15083-. 49191*1, x4 = -.56975-1.8582*1] 
[xl = .41488+. 24720*1, x2 = -1.7788+1.0599*1, 
x3 = -.41488-. 24720*1, x4 = 1.7788-1.0599*1] 
»> component #2: 

[xl = -.95775+. 36799*1, x2 = -. 90980- . 34957*1 , 
x3 = .95775-. 36799*1, x4 = . 90980+ . 34957*1] 
[xl = .71538+. 12328*1, x2 = 1 . 3576- . 23395*1 , 
x3 = -.71538-. 12328*1, x4 = -1 . 3576+ . 23395*1] 
[xl = -3.7686+1.7072*1, x2 = -. 22017- . 99738e-l*I , 
x3 = 3.7686-1.7072*1, x4 = . 22017+ . 99738e-l*I] 
[xl = -. 16036-. 30943*1, x2 = -1.3202+2.5476*1, 
x3 = .16036+. 30943*1, x4 = 1.3202-2.5476*1] 
»> component #3: 

[xl = -1.0-.53734e-17*I, x2 = 1 . 0- . 20045e-16*I , 
x3 = 1.0+.89149e-17*I, x4 = -1 . 0+ . 18026e-17*I] 

»> component #10: 

[xl = -.593516-17+1.0*1, x2 = - . 46995e-16+l . 0*1 , 
x3 = .161586-16-1.0*1, x4 = . 22439e-16-l . 0*1] 

Note that there are 2 witness sets of 4 points corresponding to the 2 isolated curves 
and 8 singletons for the embedded points. 

See the webpage [TO] for the scripts in Macaulay 2 and Maple (using PHCmaple pack- 
age [TT]) that perform prime decomposition and numerical irreducible decomposition, 
respectively, for the first deflation ideal in this example. 

6. Discussion and conclusion 

We consider this paper as one laying a theoretical foundation to the method that at 
this point works only on small examples. However, it is our believe that this technique 
would be able to solve problems unsolvable by purely symbolic methods in the future. 
The improvements are expected to be made both in the software and in the theory. 

We remark that the software in the area of numerical algebraic geometry is as young as 
the area itself. For the purposes of numerical irreducible decomposition there exist only 
two software options: PHCpack [20j that we use via PCHmaple pi] and Bertini [2]. The 
practical computation using the ideas in this paper is limited by the capabilities of these 
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software systems; we expect the implementations of numerical irreducible decomposition 
algorithms to improve. Both PHCpack and Bertini move towards the throughout par- 
allelization; as we argued in the introduction, easy parallelization is a crucial feature of 
numerical methods that distinguishes them from the symbolic ones. 

The future theoretical work should, in particular, concentrate on the construction of 
special homotopy methods to tackle deflated ideals (systems), which possess an obvious 
multihomogeneous structure: they are linear in the additional variables. Also, while the 
global algorithms such as NID have been well established, there are still no efficient local 
procedures, e.g., for determining the local dimension at a given point on a variety. The 
same can be said about the local dual space computation: while it is possible to compute 
the truncation at some degree using the deflation matrix, an efficient description of the 
whole (possibly infinite-dimensional) dual space and ways to create such description are 
yet to be found. 
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